Pegasus Tutorial¶
Compiling date: 2020-08-11
import pegasus as pg
Count Matrix File¶
For this tutorial, we provide a count matrix dataset on Human Bone Marrow with 8 donors stored in zarr format (with file extension ".zarr.zip").
You can download the data at https://storage.googleapis.com/terra-featured-workspaces/Cumulus/MantonBM_nonmix_subset.zarr.zip.
This file is achieved by aggregating gene-count matrices of the 8 10X channels using PegasusIO, and further filtering out cells with fewer than $100$ genes expressed. Please see here for how to do it interactively.
Now load the file using pegasus read_input function:
data = pg.read_input("MantonBM_nonmix_subset.zarr.zip")
data
The count matrix is managed as a UnimodalData object defined in PegasusIO module, and users can manipulate the data from top level via MultimodalData structure, which can contain multiple UnimodalData objects as members.
For this example, as show above, data is a MultimodalData object, with only one UnimodalData member of key "GRCh38-rna", which is its default UnimodalData. Any operation on data will be applied to this default UnimodalData object.
UnimodalData has a structure similar to AnnData defined in anndata package (see anndata.AnnData for details; the figure below is also borrowed from this link):
It has 5 major parts:
- Raw count matrix:
data.X, a Scipy sparse matrix, with rows the cell barcodes, columns the genes/features:
data.X
This dataset contains $48,219$ barcodes and $36,601$ genes.
- Cell barcode attributes:
data.obs, a Pandas data frame with barcode as the index. For now, there is only one attribute"Channel"referring to the donor from which the cell comes from:
data.obs.head()
data.obs['Channel'].value_counts()
- Gene attributes:
data.var, also a Pandas data frame, with gene name as the index. For now, it only has one attribute"gene_ids"referring to the unique gene ID in the experiment:
data.var.head()
- Unstructured information:
data.uns, a Python hashed dictionary. It usually stores information not restricted to barcodes or features, but about the whole dataset, such as its genome reference and modality type:
data.uns['genome']
data.uns['modality']
- Finally, embedding attributes on cell barcodes:
data.obsm; as well as on genes,data.varm. We'll see it in later sections.
pg.qc_metrics(data, percent_mito=10)
The metrics considered are:
- Number of genes: keep cells with $500 \leq \text{# Genes} < 6000$ (Default);
- Number of UMIs: don't filter cells due to UMI bounds (Default);
- Percent of Mitochondrial genes: keep cells with percent $< 10\%$.
For details on customizing your own thresholds, see documentation.
Numeric summaries on filtration on cell barcodes and genes can be achieved by get_filter_stats method:
df_qc = pg.get_filter_stats(data)
df_qc
The results is a Pandas data frame on samples.
You can also check the QC stats via plots. Below is on number of genes:
pg.qcviolin(data, plot_type='gene', dpi=100)
Then on number of UMIs:
pg.qcviolin(data, plot_type='count', dpi=100)
On number of percentage of mitochondrial genes:
pg.qcviolin(data, plot_type='mito', dpi=100)
Now filter cells based on QC metrics set in qc_metrics:
pg.filter_data(data)
You can see that $35,465$ cells ($73.55\%$) are kept.
Moreover, for genes, only those with no cell expression are removed. After that, we identify robust genes for downstream analysis:
pg.identify_robust_genes(data)
The metric is the following:
- Gene is expressed in at least $0.05\%$ of cells, i.e. among every 6000 cells, there are at least 3 cells expressing this gene.
Please see its documentation for details.
As a result, $25,653$ ($70.09\%$) genes are kept. Among them, $17,516$ are robust.
We can now view the cells of each sample after filtration:
data.obs['Channel'].value_counts()
Normalization and Logarithmic Transformation¶
After filtration, we need to first normalize the distribution of cells w.r.t. each gene to have the same count (default is $10^5$, see documentation), and then transform into logarithmic space by $log(x + 1)$ to avoid number explosion:
pg.log_norm(data)
For the downstream analysis, we may need to make a copy of the count matrix, in case of coming back to this step and redo the analysis:
data_trial = data.copy()
Highly Variable Gene Selection¶
Highly Variable Genes (HVG) are more likely to convey information discriminating different cell types and states. Thus, rather than considering all genes, people usually focus on selected HVGs for downstream analyses.
You need to set consider_batch flag to consider or not consider batch effect. At this step, set it to False:
pg.highly_variable_features(data_trial, consider_batch=False)
By default, we select 2000 HVGs using the pegasus selection method. Alternative, you can also choose the traditional method that both Seurat and SCANPY use, by setting flavor='Seurat'. See documentation for details.
We can view HVGs by ranking them from top:
data_trial.var.loc[data_trial.var['highly_variable_features']].sort_values(by='hvf_rank')
We can also view HVGs in a scatterplot:
pg.hvfplot(data_trial, dpi=200)
In this plot, each point stands for one gene. Blue points are selected to be HVGs, which account for the majority of variation of the dataset.
Principal Component Analysis¶
To reduce the dimension of data, Principal Component Analysis (PCA) is widely used. Briefly speaking, PCA transforms the data from original dimensions into a new set of Principal Components (PC) of a much smaller size. In the transformed data, dimension is reduced, while PCs still cover a majority of the variation of data. Moreover, the new dimensions (i.e. PCs) are independent with each other.
pegasus uses the following method to perform PCA:
pg.pca(data_trial, robust=True)
By default, pca uses:
- Before PCA, scale the data to standard Normal distribution $N(0, 1)$, and truncate them with max value $10$;
- Number of PCs to compute: 50;
- Apply PCA only to highly variable features.
In addition, I set robust=True for reproducibility purpose. In practice, the default setting uses randomized PCA for fast computation.
See its documentation for customization.
To explain the meaning of PCs, let's look at the first PC (denoted as $PC_1$), which covers the most of variation:
coord_pc1 = data_trial.uns['PCs'][:, 0]
coord_pc1
This is an array of 2000 elements, each of which is a coefficient corresponding to one HVG.
With the HVGs as the following:
data_trial.var.loc[data_trial.var['highly_variable_features']].index.values
$PC_1$ is computed by
\begin{equation*} PC_1 = \text{coord_pc1}[0] \cdot \text{HES4} + \text{coord_pc1}[1] \cdot \text{ISG15} + \text{coord_pc1}[2] \cdot \text{TNFRSF18} + \cdots + \text{coord_pc1}[1997] \cdot \text{RPS4Y2} + \text{coord_pc1}[1998] \cdot \text{MT-CO1} + \text{coord_pc1}[1999] \cdot \text{MT-CO3} \end{equation*}Therefore, all the 50 PCs are the linear combinations of the 2000 HVGs.
The calculated PCA count matrix is stored in the obsm field, which is the first embedding object we have
data_trial.obsm['X_pca'].shape
For each of the $35,465$ cells, its count is now w.r.t. 50 PCs, instead of 2000 HVGs.
Nearest Neighbors¶
All the downstream analysis, including clustering and visualization, needs to construct a k-Nearest-Neighbor (kNN) graph on cells. We can build such a graph using neighbors method:
pg.neighbors(data_trial)
It uses the default setting:
- For each cell, calculate its 100 nearest neighbors;
- Use PCA matrix for calculation;
- Use L2 distance as the metric;
- Use hnswlib search algorithm to calculate the approximated nearest neighbors in a really short time.
See its documentation for customization.
Below is the result:
print(f"Get {data_trial.uns['pca_knn_indices'].shape[1]} nearest neighbors (excluding itself) for each cell.")
data_trial.uns['pca_knn_indices']
data_trial.uns['pca_knn_distances']
Each row corresponds to one cell, listing its neighbors (not including itself) from nearest to farthest. data_trial.uns['pca_knn_indices'] stores their indices, and data_trial.uns['pca_knn_distances'] stores distances.
Clustering and Visualization¶
Now we are ready to cluster the data for cell type detection. pegasus provides 4 clustering algorithms to use:
louvain: Louvain algorithm, using louvain package.leiden: Leiden algorithm, using leidenalg package.spectral_louvain: Spectral Louvain algorithm, which requires Diffusion Map.spectral_leiden: Spectral Leiden algorithm, which requires Diffusion Map.
See this documentation for details.
In this tutorial, we use the Louvain algorithm:
pg.louvain(data_trial)
As a result, Louvain algorithm finds 19 clusters:
data_trial.obs['louvain_labels'].value_counts()
We can check each cluster's composition regarding donors via a composition plot:
pg.compo_plot(data_trial, 'louvain_labels', 'Channel')
However, we can see a clear batch effect in the plot: e.g. Cluster 10 and 13 have most cells from Donor 3.
We can see it more clearly in its FIt-SNE plot (a visualization algorithm which we will talk about later):
pg.fitsne(data_trial)
pg.scatter(data_trial, attrs=['louvain_labels', 'Channel'], basis='fitsne')
Batch Correction¶
Batch effect occurs when data samples are generated in different conditions, such as date, weather, lab setting, equipment, etc. Unless informed that all the samples were generated under the similar condition, people may suspect presumable batch effects if they see a visualization graph with samples kind-of isolated from each other.
For this dataset, we need the batch correction step to reduce such a batch effect, which is observed in the plot above.
In this tutorial, we use Harmony algorithm for batch correction. It requires redo HVG selection, calculate new PCA coordinates, and apply the correction:
pg.highly_variable_features(data, consider_batch=True)
pg.pca(data, robust=True)
pca_key = pg.run_harmony(data)
The corrected PCA coordinates are stored in data.obsm:
data.obsm[f"X_{pca_key}"].shape
Repeat Previous Steps on the Corrected Data¶
As the count matrix is changed by batch correction, we need to recalculate nearest neighbors and perform clustering. Don't forget to use the corrected PCA coordinates as the representation:
pg.neighbors(data, rep=pca_key)
pg.louvain(data, rep=pca_key)
Let's check the composition plot now:
pg.compo_plot(data, 'louvain_labels', 'Channel')
If everything goes properly, you should be able to see that no cluster has a dominant donor cells. Also notice that Louvain algorithm on the corrected data finds 16 clusters, instead of the original 19 ones.
Also, FIt-SNE plot is different:
pg.fitsne(data, rep=pca_key)
pg.scatter(data, attrs=['louvain_labels', 'Channel'], basis='fitsne')
You can see that the right-hand-side plot has a much better mixture of cells from different donors.
Visualization¶
tSNE Plot¶
In previous sections, we have seen data visualization using FIt-SNE. FIt-SNE is a fast implementation on tSNE algorithm. Including FIt-SNE, pegasus provides 3 different tSNE plotting methods:
fitsne: FIt-SNE plot, using fitsne package. See detailstsne: tSNE plot, using Multicore-TSNE package. See detailsnet_tsne: approximated tSNE plot with speed up using Deep Neural Networks (DNN) models. See details
UMAP Plot¶
Besides tSNE, pegasus also provides UMAP plotting methods:
umap: UMAP plot, using umap-learn package. See detailsnet_umap: Approximated UMAP plot with DNN model based speed up. See details
Below is the UMAP plot of the data using umap method:
pg.umap(data, rep=pca_key)
pg.scatter(data, attrs=['louvain_labels', 'Channel'], basis='umap')
Differential Expression Analysis¶
With the clusters ready, we can now perform Differential Expression (DE) Analysis. DE analysis is to discover cluster-specific marker genes. For each cluster, it compares cells within the cluster with all the others, then finds genes significantly highly expressed (up-regulated) and lowly expressed (down-regulated) for the cluster.
Now use de_analysis method to run DE analysis. We use Louvain result here.
pg.de_analysis(data, cluster='louvain_labels', auc=False)
By default, DE analysis runs
auc: Area under ROC (AUROC) and Area under Precision-Recall (AUPR).t: Welch’s t test.
Alternatively, you can also run the follow tests by setting their corresponding parameters to be True:
fisher: Fisher’s exact test.mwu: Mann-Whitney U test.
DE analysis result is stored with key "de_res" (by default) in varm field of data. See documentation for more details.
To load the result in a human-readable format, use markers method:
marker_dict = pg.markers(data)
By default, markers:
- Sort genes by WAD scores in descending order;
- Use $\alpha = 0.05$ significance level on q-values for inference.
See documentation for customizing these parameters.
Let's see the up-regulated genes for Cluster 1, and rank them in descending order with respect to log fold change:
marker_dict['1']['up'].sort_values(by='log_fold_change', ascending=False)
Among them, TRAC worth notification. It is a critical marker for T cells.
We can also use Volcano plot to see the DE result. Below is such a plot w.r.t. Cluster 1 with log fold change being the metric:
pg.volcano(data, cluster_id = '1', dpi=200)
The plot above uses the default thresholds: log fold change at $1$ (i.e. fold change at $2$), and q-value at $0.05$. Each point stands for a gene. Red ones are significant marker genes: those at right-hand side are up-regulated genes for Cluster 1, while those at left-hand side are down-regulated genes.
We can see that gene TRAC is the second to rightmost point, which is a significant up-regulated gene for Cluster 1.
To store a specific DE analysis result to file, you can write_results_to_excel methods in pegasus:
pg.write_results_to_excel(marker_dict, "MantonBM_subset.de.xlsx")
Cell Type Annotation¶
After done with DE analysis, we can use the test result to annotate the clusters.
celltype_dict = pg.infer_cell_types(data, markers = 'human_immune', de_test = 't')
cluster_names = pg.infer_cluster_names(celltype_dict)
infer_cell_types has 2 critical parameters to set:
markers: Either'human_immune','mouse_immune','human_brain','mouse_brain','human_lung', or a user-specified marker dictionary.de_test: Decide which DE analysis test to be used for cell type inference. It can be either't','fisher', or'mwu'.
infer_cluster_names by default uses threshold = 0.5 to filter out candidate cell types of scores lower than 0.5.
See documentation for details.
Below is the cell type annotation report for Cluster 1:
celltype_dict['1']
The report has a list of predicted cell types along with their scores and support genes for users to decide.
Next, substitute the inferred cluster names in data:
- Construct an annotation dictionary of the following format:
anno_dict = dict(zip(map(lambda n: str(n), range(1, len(cluster_names)+1)), cluster_names))
anno_dict
The annotation dictionary has keys being cluster labels, and their putative cell types being values. In practice, users may want to manually create this dictionary by reading the report in celltype_dict.
- Use
annotatefunction to rename cluster labels:
pg.annotate(data, name='anno', based_on='louvain_labels', anno_dict=anno_dict)
data.obs['anno'].value_counts()
I use "anno" for the key name of the new cluster names.
Now plot data with cluster names:
pg.scatter(data, attrs='anno', basis='fitsne', dpi=100)
pg.scatter(data, attrs='anno', basis='umap', legend_loc='on data', dpi=150)
Raw Count vs Log-norm Count¶
Now let's check the count matrix:
data
You can see that besides X, there is another matrix raw.X generated for this analysis. As the key name indicates, raw.X stores the raw count matrix, which is the one after loading from the original Zarr file; while X stores the log-normalized counts.
data currently binds to matrix X. To use the raw count instead, type:
data.select_matrix('raw.X')
data
Now data binds to raw counts.
Gene-specific Plots¶
You may want to check expression of a specific group of genes. Pegasus provides a number of plots on this purpose.
Let's check the following genes for example:
marker_genes = ['CD38', 'JCHAIN', 'FCGR3A', 'HLA-DPA1', 'CD14', 'CD79A', 'MS4A1', 'CD34', 'TRAC', 'CD3D', 'CD8A',
'CD8B', 'GYPA', 'NKG7', 'CD4', 'SELL', 'CCR7']
Also, since most of the plots below use log-norm counts, I'll select matrix X for convenience:
data.select_matrix('X')
Violin Plot¶
pg.violin(data, attrs=marker_genes, groupby='anno', dpi=100)
Heatmap¶
By default, Pegasus uses average expressions within clusters for heatmap plotting:
pg.heatmap(data, genes=marker_genes, groupby='anno', row_cluster=True)
You can also show expressions of individual cells:
pg.heatmap(data, genes=marker_genes, groupby='anno', on_average=False)
Dotplot¶
pg.dotplot(data, genes=marker_genes, groupby='anno')
Feature Plot¶
Feature plot uses the same method scatter when plotting cell embeddings, except that this time attrs are fed with gene names.
pg.scatter(data, attrs=['TRAC', 'CD79A', 'CD14', 'CD34'], basis='umap')
Dendrogram¶
Below shows hierarchical clustering result using expressions regarding marker_genes on clusters:
pg.dendrogram(data, genes=marker_genes, groupby='anno', dpi=100)
You can also plot dendrogram based on the corrected PCA coordinates of cells:
pg.dendrogram(data, groupby='anno', rep=pca_key, dpi=100)
Please see Plotting documentation for details on these plotting methods.
Save Result to File¶
Use write_output function to save analysis result data to file:
pg.write_output(data, "result.zarr.zip")
It's stored in zarr format, because this is the default file format in Pegasus.
Alternatively, you can also save it in h5ad, mtx, or loom format. See its documentation for instructions.